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The phase space density of fermionic dark matter haloes 
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ABSTRACT 

We have performed a series of numerical experiments to investigate how the primordial ther- 
mal velocities of fermionic dark matter particles affect the physical and phase space density 
profiles of the dark matter haloes into which they collect. The initial particle velocities in- 
duce central cores in both profiles, which can be understood in the framework of phase space 
density theory. We find that the maximum coarse-grained phase space density of the simu- 
lated haloes (computed in 6 dimensional phase space using the EnBid code) is very close to 
the theoretical fine-grained upper bound, while the pseudo phase space density, Q ~ p/cr 3 , 
overestimates the maximum phase space density by up to an order of magnitude. The density 
in the inner regions of the simulated haloes is well described by a 'pseudo-isothermal' pro- 
file with a core. We have developed a simple model based on this profile which, given the 
observed surface brightness profile of a galaxy and its central velocity dispersion, accurately 
predicts its central phase space density. Applying this model to the dwarf spheroidal satellites 
of the Milky Way yields values close to 0.5 keV for the mass of a hypothetical thermal warm 
dark matter particle, assuming the satellite haloes have cores produced by warm dark matter 
free streaming. Such a small value is in conflict with the lower limit of 1.2 keV set by obser- 
vations of the Lyman-a forest. Thus, if the Milky Way dwarf spheroidal satellites have cores, 
these are likely due to baryonic processes associated with the forming galaxy, perhaps of the 
kind proposed by Navarro, Eke and Frenk and seen in recent simulations of galaxy formation 
in the cold dark matter model. 
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5t ■ 1 INTRODUCTION 

The standard cosmological model, the "A cold dark matter model" 
(ACDM) has been tested over a huge range of scales, from the 
entire observable universe, probed by measurements of tempera- 
ture an isotropies in the cos mic microwave background radiation 
(CMB; [Komatsu et alj|201 lb . to the scales of galaxy clusters and 
individual bright galaxies, probed by large galaxy and Lyman-q 
forest surveys dColless et al. 2005; Se liak et alj|2005j : IZehavi et al.l 
l25Tlh . On smaller scales than this, there is no strong evidence to 
support the standard model. Yet, it is on such scales that the nature 
of the dark matter is most clearly manifest. In the standard model 
the dark matter consists of cold particles, such as the lightest stable 
particle predicted by Supersymmetry. There are, however, models 
of particle physics that predict lighter particles, such as sterile neu- 
trinos, that would behave as warm (WDM), rather than cold dark 
matter (CDM) (see Feng 2oT(il: lHoopeill2012l for discussions of re- 
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cent experimental constraints). No current astronomical data can 
distinguish between these alternatives. 

If the dark matter particles freeze out while in thermal equi- 
librium, their kinematics leave a pred ictable imprint in the power 
spectrum of primordial perturbations t Zerdovichll 1965 ). Hot dark 
matter particles, such as light neutrinos, decouple while still rela- 
tivistic; their large thermal velocities dampen fluctuations below a 
'free streaming' length, Afs, which is of the order of a f ew tens 
of megaparsecs at redshift z — dBond & Szalavl lT983). Early 
N-body simulations showed that this is too large to be compati- 
ble with the level of c lustering measured in the galaxy distribu- 
tion dFrenketal.ll 19831) , thus ruling out light neutrinos as the dom- 
inant form of dark matter. Recent analyses o f the power spectrum 
of the galaxy di stribution dCole et aL 20051). the Lv man-q forest 
dViel et al.ll2010h and the CMB dKomatsu et alj|201ll) constrain the 
sum of the neutrino masses to be ^ rn v < 0.58 eV. 

Cold dark matter particles decouple after they have become 
non-relativistic. For a typical cold, weakly interacting, massive par- 
ticle, Afs is of the order of a par sec and the corre sponding 'Jeans' 
mass is of the order of 1O" 6 M dGreen et alj2 005); free streaming 
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in this case is not relevant for galaxy formation. The intermediate 
case of warm dark matter corresponds to particles that decouple 
while still relativistic, yet become non-relativistic before the epoch 
of radiation-matter equality. Their free streaming scale could then 
be of order the size of a galaxy, in which case they can affect the 
build- up of dark matte r haloes and of the galaxies forming within 



them jBode et alJ200ll;lLovell et al.l2012l) . as well as the formation 
of the first stars IGao & Theunsl 2007 ). Dark matter particles need 
not freeze out in thermal equilibrium, as is the ca se, for example, 
of the a xion and the sterile neutrino discussed by IBovarskv et al.l 
(2009a). In this case, the relation between the dark matter particle 
mass and Afs is more complex. IBovarskv etall (2009b) discuss 
current limits on such WDM models. 

Except perhaps in the central regions of galaxies, the dynam- 
ics of dark matter particles are driven purely by their own gravity 
and are thus governed by the Vlasov-Poisson equations. Then, ac- 
cording to Liouville's theorem, the distribution function, /(x, v, t), 
of the particles - the mass per unit volume in phase space - is time- 
independent, -D/(x, v,t)/ Dt — 0. Phase space mixing, or coarse- 
graining, can only decrease the phase space densi ty below this fine- 
grain bound set by the nature of the particles. iTremaine & Gunr] 
( 1979) used this property to constrain the nature of stable leptons 
as the dominant form of dark matter. Since the intrinsic velocities 
of CDM particles are small, their fine-grained phase s pace density 
limit is very high and such haloes are cuspy (e.g. 
1996tll997l:lDiemand et al. 2008: IStadel et alj|2009l : 
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20101 : iGao et al.ll2012h . For WDM (and a fortiori for HDM), the 
phase space density bound is much lower and such haloes develop 
central cores (Hog an&Dalcanton 2000; Bo de et alj200lT) . 

Evidence for ce ntral cores in galaxies is controversial (see 
IFrenk & White] l20ll and references therein for a recent dis- 
cussion). Recently, Walker & Penarrubial d2012h have argued that 
the kinematics of the Fornax and Sculptor dwarf sp heroidals in 
the Mi lky Way rule out central density cusps, but Istrigari et al.1 
(2010a) have shown that the data for these and other dwarfs 
are consistent with the cuspy pro files seen in cold da r k mat - 
ter simulations. On the other ha nd, Bovla n-Kolchin et "all feOllI) . 
ILovell et all d2012l) and IParrv et alj J2012h have shown that the 
central concentration of the dark matter haloes of nine dwarf 
spheroidals ar e lower than expecte d in sim ulations of cold da rk 
matter haloes dSpringel et afl fcoOS) (but see IWang et all d2012l) ). 
Although baryonic processes associated with the forming galaxy 
could lower the cent r al density of haloes and even induce a core 
dNavarro et al.1 1 19961; feead & Gilmord 120051 ; iMashchenko et al.1 
120081 : iPontzen & Gove rnato 2012), it is also possible, in principle, 
that the kind of phase space constraints just discussed could be at 
work at the centres of the dwarf spheroidals. 

In this paper we investigate the effect of primordial dark mat- 
ter particle thermal velocities on halo density profiles. To this aim, 
we carry out N-body simulations from initial conditions in which 
the power spectrum has a small-scale cut-off and the particles rep- 
resenting the dark matter have significant thermal velocities. Our 
setup does not correspond exactly to any particular HDM or WDM 
particle candidate. Rather, we are interested in the more general 
problem of how phase space constraints are satisfied in cosmo- 
logical N-b ody simulations. In particu lar, we determine the level 
at which the Tremain e & Gunnl Jl979i) bound is satisfied. We then 
model our numerical results, generalise them to the case of specific 
WDM candidates, and apply them to Milky Wa y satellites. 

As we were completing this work. lMaccio et ai]d2012l) posted 
a paper on Arxiv investigating similar issues to those of interest 
here, using similar numerical techniques. However, a numerical er- 



ror in their initial condition^] requires a re-interpretation of their 
results. 

The remainder of this paper is structured as follows. In Sec- 
tion 2, we describe our set of numerical simulations. In Section 3, 
we briefly review relevant aspects of phase space theory and de- 
scribe our methods for calculating the coarse-grained phase space 
density in the simulations. In Section 4, we investigate phase and 
real space density profiles of dark matter haloes and compare these 
with our theoretical estimates. In Section 5, we use a model based 
on our numerical simulations to revise the lower mass limit of 
WDM particles. Our paper concludes in Section 6 with a summary 
and discussion. 



2 PHASE SPACE DENSITY 

In this section, we briefly review those aspects of phase space den- 
sity theory that are relevant to general fermionic dark matter parti- 
cles. 



2.1 Fine-grained phase space density 

The evolution of a system of collisionless particles is described by 
the Vlasov equation. According to Liouville's theorem, the fine- 
grained phase space density, /(x, v, t) - the mass density in an in- 
finitesimal six-dimensional phase space volume, d 3 xd 3 v, centred 
on the point (x, v) a t time t - is conserved: Df/Dt = (see e.g. 

iPeebles & Yul d 19701) for the General R elativi stic descripti on). 

Following ITremaine & Gunnl dl979l) . iMadserJ dl99ll) 
an dHogan & Dalcantonl d2000h , using the Fermi-Dirac occupation 
distribution we can derive an upper limit on J-fd, the fine-grained 
phase space density of a relativistic fermionic relic in kinetic 
equilibrium: 



Jfd(p) 



9 

2(2nh) 3 



,(B- P )/(fc B r) _|_ i 



(1) 



Here, E = [(rac 2 ) 3 + (pc) 2 ] 1 ^ 2 is the energy of a particle 
with mass m and momentum p, and the chemical potential, 
fi = for thermal relics ; g is the number of degrees of free- 
dom dKolb & TurneTl 990). This expression may be more familiar 
from the derivation of the equation of state of degenerate fermions. 
Warm dark matter particles are relativistic when they decouple, i.e. 
when T = To, the decoupling temperature. After decoupling both 
the energy and the temperature decline in proportion to (1 + z) 
and hence the shape of the Fermi-Dirac distribution remains un- 
changed. Finally, under a Lorentz transformation, lengths contract 
in proportion to the Lorentz factor, 7, whereas momenta increase 
in proportion to I/7 so the phase space density is an invariant. 

The function Jfd is the number density of particles per unit 
volume in momentum space, an appropriate choice when the parti- 
cles are relativistic. When the particles become non-relativistic, we 
can write the phase space density, /fzj, in terms of the mass den- 
sity of particles per unit volume in velocity space, which introduces 



1 The velocity dispersion o f the dark matter particles in the initial condi- 
tions of Maccio et al. (2012) is a factor of \/3 larger than assumed. 
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a factor m 4 , so the non-relativistic version of Eqn. (TJ becomes 

4 

rmax Q 

= OA2M ^- 3 0™*-r 3 l(^f- (2) 

In the linear regime, the density of these non-relativistic particles 
decays as (1 + z) 3 but their peculiar velocities decay as (1 + z); 
hence this bound does not evolve with redshift, as expected. 



To derive the limit in Eqn. l[5) we had to make an assump- 
tion about the form o f the density profile of the central halo. 
iBovarskv et al.l d2009ah introduced a 'coarsest' phase space den- 
sity, FgJ,y X , which corresponds to a maximal coarse graining and 
does not require this assumption. Instead it makes use of the value 
of the enclosed mass, M(R), of the halo over the whole available 
phase space volume Ax Av = (47r/3) 2 i? 3 «? M , where is the 
local escape speed. For isotropic orbits, Voo > v6u, leading to the 
maximum coarse-grained value of F of 



2.2 Coarse-grained phase space density 

The coarse-grained phase space density, F(pc,\,t), is defined as 
the mass density in a. finite six-dimensional phase space volume, 
A 3 xA 3 v, centered on the point (x, v) at time t. The averaging 
process decr eases the phase space densi ty and hence F < f (see, 
for example, iTremaine et ail dl986h and lMathun dl988l) for the ap- 
plication of Liouville's theorem in an astronomical setting). Phase 
mixing is expected to be small in the region with the highest phase 
space density, and the maximum value of F for a system is th ere- 
fore expected to be close to the value of / dLvnden-B ellll 1 96% . Of 
course, in any real system, baryonic effects may increase or de- 
crease the phase space density. In a simulated dark matter halo 
we can estimate F in 6 dimensions and verify whether indeed 
F m " < fiffi*; we will do so below. 

We begin by discussing various ways in which F can be es- 
timated for a dark matter halo. We assume the centra l part of the 
halo to have a pseud o-isothermal density profile (e.g. iKentll 19861 : 
iBegeman et al.ll99lh , which is a good fit to the simulated dark mat- 
ter halo discussed below. This profile is given by 



p(r) = 



Po 



1 + {r/r c y 



(3) 



where po ' s the core density and r c the core radius. The correspond- 
ing asymptotically flat circular velocity is V c = (47rGpo?"c) 1 ' /2 - 
Assuming isotropic orbits V c is related to the one-dimensional ve- 
locity dispersion a, by V c = V2a. Under these assumptions, cen- 
tral density, core radius and velocity dispersion are related by 



po 



2-kG r c 



(4) 



To obtain the corresponding coarse-grained phase space den- 
sity, we need to know the velocity distribution function. Dynam- 
ically relaxed systems often have a Maxwellian velocity distribu- 
tion, and we will assume this to be a good description for the cen- 
tral regions of the dark matter halo ( Vogelsberger et al. 2009). The 
maximum density in velocity space is then (27rcr 2 ) _3//2 , and there- 
fore the maximum coarse-grained phase space density, which oc- 
curs at the centre of the halo, is 
Po 



(2na 2 ) 3 / 2 
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(2tt) 5 / 2 G^ 
0.06 Mq kpc" 3 (kms -1 ) 



x (- 



(5) 



(6) 



100 km s" 1 ' '20 kpc' 

If the core is due to free-streaming of a WDM particle, then requir- 
ing F™ ax < /™d x constrains the WDM particle mass to be 



mc 2 (f )^ 4 > 8.2 keV (^-^y 1 ^ (^)" 1/2 



pc 



(7) 
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0.18F™ ax , (8) 



87rV6cr 3 16^tt 2 G or\ 

where is the half-light radius, defined by IPrvor & Kormendvl 
dl990h . 



Po 



3 In 2 a 1 
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(9) 



This calculation is less restrictive and hence -FgJ, ax < FS?*< as 
should be. 

Finally, an often used approximation to F is 



Q 



(10) 



introduced by iHogan & Dalcantonl j2000j)_ and more recently 



used, fo r example, by [ Taylor & Navarrol i 200 lh: Ascasibar et al.l 
d2004h; iDe hnen & McLaughlin j2005h; IPeirani et alj d2006t); 



[ Hoffman et all d2007l) : IColmetalJ d2008h ; IVass et all d2009h : 
iNavarro et al.l d2010h . We will refer to this as a pseudo phase space 
density, since it only has the dimensions of F, but it is not a proper 
coarse-grained average of /. Indeed, the maximum value of Q 
for the halo profile of Eqn. l|3j under the assumptions of isotropic 
Maxwellian velocities, is 



Q max = 



(27T) 



3/2 



3^3 



-7^ ax «3.03F£r 



(11) 



The upper bound to the fine-grained dis t ributio n for thermal 
fermions discussed by Hoga n & Dalcantonl d200d) is very simi- 
lar to the value in Eqn .Q above, /noTan ~ 0.97/™d x (see also 
IBovarskv et al. I d2009j) ). Requiring Q max < /^ ax an then overesti- 
mates the constraint on the dark matter particle mass by a factor of 
(3.12) 1/4 . 



2.3 Coarse-grained phase space density in simulations 

A robust measurement of F in an iV-body system requires 
a full six-dimensional calculation. In recent years, a num- 
ber of independent a pproaches have been developed that allow 
such a calculation Arad et al] |2004|; lAscasibar & Binnevl 120051 : 
ISharma & Steinmetzl l2006; Vogelsber ger et alj|2008l) . The method 
of lAradetal.ld2004l) uses the Delaunay tessellation of the N parti- 
cles in 6 dimensional (x, v) phas e space. Tw o other algorithms, Fi- 
EstAS d Ascasibar & Binney|2 005) and EnBiD dSharma & Steinmetj 
2006), are based upon the idea of a binary k — d-tree, i.e. repeated 
subdivisions of phase space in each of its 6 dimensions into nodes 
that contain (approximately) the same number of particles until 
each node contains a single particle. Both the Delaunay tessellation 
and the k — d tree then associate a 'volume,' Vi , to a particle, and the 
corresponding phase space density is simply m P artici c /V5, where 
wiparticic is the mass of the particle in the simulation. These algo- 
rithms contain adjustable pa rameters that can be chosen to improve 
the estimates. For example, lAscasibar & Binnevl d2005h include a 
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Figure 1. Projected density in a 5h~ x Mpc cube at z = for the CDM, WDM, and P-WDM simulations (left to right). The top panels show the full simulation 
volume in a slice of 1.5h -1 Mpc ; the bottom panels zoom in on the most massive halo. As expected, features on scales larger than the free-streaming length 
are very similar, but on smaller scales the WDM runs have much less substructure. Model P-WMD, which had substantial initial random velocities, has almost 
no substructure at all. 



boundary correction to improve the behaviour of the method at low 
density and smoothing kernels to reduce particle noise. 

I n this paper we make use of the publicly available EnBiD 
code 1 Shar ma~& Steinmetzl l200rj) to estimate the coarse-grained 
phase space density in our simulations. EnBiD employs a Shannon 
Entropy formulation which gives accurate results in high density 
regions. We have included the 'boundary correction' but did not 
use the smoothing option because it may underestimate the highest 
phase space density in a dark matter halo, which is crucial in this 
study. We will denote the coarse-grained phase space density of the 
simulated haloes calculated using EnBiD as Fnb 



3 iV-BODY SIMULATIONS 
3.1 Physical requirements 

If dark matter particles have large intrinsic velocities at early times, 
two distinct physical effects become important: free streaming out 
of density perturbations and a maximum achievable phase space 
density. In the linear regime, the WDM power spectrum of fluctu- 
ations, Pwdm(&), is suppressed relative to the CDM power spec- 
trum, PcDM(fe), by a factor 



T 2 (fc) 



tW 



-PCDM (fc) 



(12) 



where k is the wavevector. Free streaming introduces a smallest 
wavevector, /ch, for which T < 1, and some authors refer to 2n/ku 
as the 'free-streaming length'. However, when considering struc- 
ture formation, it is more useful to characterise the effects of free 
streaming by the value of ftps = 27t/Afs for which the amplitude 
of the power spectrum is reduced by a factor two relative to the 
CDM case, such that T 2 (k = k FS ) = 1/2. The shape of the func- 
tion T(k) depends on the nature of the dark matter. A commonly 
used approximation for thermally produced WDM particles is 



T(k) = (l + (ak) 2 Y 
jBode et alfcOOll) . For this case, 



(13) 



A FS = 23.45 a, 



-fiv 



r iWDM .0.15/ n sl.3/ mc \ 



'/i _1 Mpc, (14) 



where Owdm is the total WDM density today in units of 
the critical density, m is the WDM particle mass, and H = 
100 h km s _1 Mpc -1 is the Hubble constant at z — (Bo de et all 
1200 lb . The numerical coefficient depends on the effective number 
of degrees of freedom, gx, of the particle and we have assumed 
gx = 1.5. With this definition we find 



A FS = 1.2 ( 



f2wDM No.iSf h ,i.3, mc . 



0.3 



'keV 



5 /i _1 Mpc; (15) 
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Figure 2. Mean, radially averaged density profile, p(r) (left), and coarse-grained phase space density profile, -FnbM (right), for the most massive halo in 
the CDM (blue), WDM (green) and P-WDM (red) simulations. Phase space densities were computed using the EnBiD code. The black-dashed lines are NFW 
and power-law fits to p and -Fnb respectively for the CDM model. Different line styles correspond to simulations of different numerical resolution, as given 
in the legend, with arrows indicating twice the corresponding softening lengths. The CDM and WDM profiles are nearly indistinguishable, with p(r) and 
FNg(r) well fit by the NFW and power-law profiles respectively. In contrast, the P-WDM model, which has significant initial thermal velocities, develops a 
core. The maximum value of the phase space density, -Fnb, in the P-WDM model is close to, but smaller than the fine-grained upper limit, /p^J, indicated by 
the horizontal red line. 



the corresponding mass is 

4-7T 3 

Mfs = — Pm Aps 

= 6 x 10 11 fc -i M 
v 0.3 ; v 0.72 ; v keV ; ° 

(16) 

We will refer to Afs as the free-streaming scal^l 

The second important effect is the limit on the phase space 
density achievable by dark matter particles with significant pri- 
mordial thermal velocities. This induces a core radius, r c , in col- 
lapsed haloes which, for thermal dark matter relics, scales with 
the velocity dispersion o f the halo, cr, as r c oc a~ 1 ^ 2 rn~ 2 
faogan & Dalcantonl feoOO ; Bode et al. 2001). Free streaming pre- 
vents haloes of virial radiufl R200 ~ Afs from forming in sig- 
nificant numbers. As we will show below, the core radius of such 
haloes is r c <C Afs. Since r c oc <j -1 / 2 , the ratio r c /i?2oo is even 
smaller for more massive haloes. 

The small value of the ratio r c / makes it challenging to 
investigate core radii in simulations since the calculation needs to 
resolve the huge dynamic range between r c and Afs. This is not 
currently practical for values of Afs typical of those expected in 
WDM models. However, if one is interested in the more general 

2 Unfortunately, not all authors agree on a single definition of Aps . 

3 We define the virial radius, -R200 , of a cosmological dark matter halo as 
the radius within which the mean density is 200 times the critical density. 
The corresponding enclosed mass is M200 ■ 



problem of the relationship between r c and the phase space den- 
sity, it is not necessary for the cut-off in the power spectrum to 
correspond to the thermal velocities of a particular dark matter can- 
didate. This is the approach we take in this study: we perform a nu- 
merical experiment by choosing a value of Afs but giving the par- 
ticles much larger intrinsic velocities than would be appropriate for 
a WDM particle of mass m according to Eqn. 114t . For brevity, we 
will refer to these models as pseudo-WDM models (P-WDM). We 
stress that this is not a self-consistent dark matter model: the result- 
ing core radii will be much larger than expected for a WDM particle 
of that mass, but this is imm aterial for our purpos es. This procedure 
is similar to that adopted bv lMaccio et al.l J2012h . For comparison, 
we also investigate a self-consistent WDM model in which the ini- 
tial thermal velocities are compatible with the cut-off in the power 
spectrum. Since the thermal velocities are now much smaller than 
the initial linear gravitational velocities, we expect their effect to be 
negligible. 



3.2 Simulation details 

All our simulations are of periodic cubic volumes of lin- 
ear size 5/i _1 Mpc. The assumed cosmological parameters are 
(Q m ,QA,h, as) = (0.3, 0.7, 0.7, 0.9), but the exact choice of 
these numbers is not important here. We have run two refer- 
ence simulations: a standard ACMD (labelled CDM) and a self- 
consistent A WDM (labelled WDM) model. These have 256 3 
particles, corresponding to an A-body particle mass of 6.2 x 
lO 5 h _1 M0. The CDM transfer function was computed using 
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Table 1. Simulation details. In all cases the assumed cosmological param- 
eters are (Q m ,Q^, h, ag) = (0.3,0.7,0.7,0.9) and the simulation cube 
is 5/i _1 Mpc on a side. The WDM and pseudo-WDM (P-WDM) models 
use the CDM power spectrum suppressed by T 2 (k) from Eqn. U31 . using 
the parameters for a thermal WDM particle of mass mc 2 = 1 keV with 
g = 2 spin degrees of freedom. The corresponding free-streaming length 
is Afs ~ 0.5 /i — 1 Mpc. The added velocities for the WDM model are 
consistent with the WDM particle properties, whereas for the P-WDM run, 
they correspond to those of an mc 2 = 0.03 keV thermal WDM particle. 



Name 


Number of particles 


particle mass 


softening e 






M 


kpc 


CDM 


256 3 


6.2 x 10 5 


0.8 


WDM 


256 3 


6.2 x 10 5 


0.8 


P-WDM 256 


256 3 


6.2 x 10 5 


0.8 


P-WDM512 


512 3 


7.8 x 10 4 


0.2 



CMBfast dSeljak & Zaldarriagj[l99f3) . For the WDM simulation, 
the power spectrum of the initial conditions was obtained by mul- 
tiplying the CDM spectrum by T(k) 2 , as given by Eqn. J 1 3b . The 
assumed free streaming length scale, Afs ~ 0.5 h _1 Mpc , cor- 
responds to a 2 keV thermal relic WDM particle and each particle 
is given the appropriate random velocity drawn from a Fermi-Dirac 
distribution. For definiteness we will always assume the WDM par- 
ticle has g = 2 spin degrees of freedom. 

The pseudo-WDM simulations (labelled P-WDM) have a sim- 
ilar set up to the WDM case, except that the random velocities 
are much larger and would correspond to those appropriate to a 
WDM particle of mass 0.03 keV. For ease of comparison, the Gaus- 
sian random field used to initialise the simulations uses the same 
phases for all three models. To enable statistical comparisons, we 
carried out a further 3 P-WDM simulations using the same cos- 
mological parameters but different realisations of the initial con- 
ditions; these are denoted P-WDM256- Finally, to investigate nu- 
merical convergence we carried out one simulation with 8 times 
better mass resolution, also with several different realisations of 
the initial conditions. These are denoted by P-WDM512. The start- 
ing redshift for the CDM and WDM run was z = 100, but for the 
P-WDM runs it was z = 20 since in this case structure forma- 
tion is significantly delayed compared to the CDM case. Table Q] 
summarises the parameters of the simulations. All our simulations 
were perfor med with the G ADGET- 3 code, an improved version of 
GADGET-2 dSpringelll2005l) . 

3.3 Comparisons of CDM and WDM simulations 

Images of the CDM, WDM and P-WDM256 simulations are dis- 
played in Fig. Q] for the fiducial 256 3 particle run. As anticipated, 
the overall appearance of the main halo in the three simulations 
is very similar, but the amount of small-scale substructure is very 
different. The cut-off in the initial power spectrum in the WDM 
simulations has the effect of suppressing the formation of structure 
on scales below Afs • This is the reason why there are far fewer sub- 
structures in the WDM than in the CDM images. The large initial 
random velocities in the P-WDM case further smooths out the mass 
distribution, suppressing the formation of structures with velocity 
dispersion smaller than the amplitude of the random velocities. 

The structure of the most massive halo in the CDM model, 
and the corresponding haloes in WDM and P-WDM, are compared 
in more detail in Fig. [2] These haloes all have fairly similar to- 
tal mass, M 20 o ~ 1.5 x lO 12 rT 1 M for both CDM and WDM, 



and M200 ~ 1.0 x 10 12 h _1 M Q for P-WDM512. The density pro- 
files of the CDM and WDM runs are nearly indistinguishable from 
each other, demonstrating that for this self-consistent WDM model 
neither free streaming nor the (relatively small) initial intrinsic ve- 
locities affect the overall structure of the halo. For this choice of 
WDM par ameters this is expected and has be en seen in previous 
work (e.g. lBode et~al]|200ll ; IColm et alj|2008h . In contrast, the P- 
WDM256 model, with its significantly larger intrinsic velocities, 
does form a core. Comparison of models P-WDM256 and the eight 
times higher resolution P-WDM512 shows that the core is numeri- 
cally well resolved and converged. 

The phase space density profiles of models CDM and WDM 
are also virtually indistinguishable; they are consistent with a 
power law, as found in previous studies (Fig. [2] right panel). Once 
again, neither free streaming nor intrinsic velocities affect the pro- 
file. On the other hand, in the P-WDM models a well-resolved 
core in the phase space profile develops. The coarse-grained phase 
space density, Fnb, computed using EnBiD, has a maximum of 
^nb x ~ 0.5/i 2 M & kpc -3 (km s)" 3 , close to but smaller than 
the value of = 0.85ft 2 M kpc -3 (km s) -3 (shown as the 

horizontal red line) calculated from Eqn. ((2} for a WDM particle 
of mass mc 2 — 0.03 keV. This demonstrates that our simulation 
methods and our calculation of phase space densities using EnBiD 
are sufficiently accurate to follow correctly the formation of cores 
due to the WDM intrinsic velocities. Given the good numerical con- 
vergence of the profiles, we use our three P-WDM256 simulations 
with different initial random phases to obtain statistical results in 
more detail below. 

In contrast, the fine-grained phase space density bound for the 
standard 2keV thermal relic is / F 7(mc a = 2keV) = 1.7 x 
10 7 h 2 MQ kpc -3 (kms -1 ) -3 , ~ 7 orders of magnitude larger 
than for the P-WDM case. The corresponding size of the core is 
nearly two orders of magnitude smaller, and is well within the 
smoothing length of our simulations. Clearly, as we had antici- 
pated, these simulations are not able to resolve the core of the 
WDM model. 

3.4 Physical and phase space density profiles 

As we have shown in Fig. [2] the most massive halo in the P- 
WDM simulations develops a central, near uniform density core 
of radius a few kiloparsecs. This behaviour is very different from 
the c uspy, p(r) oc 1/r, profil es familiar from CDM simula- 
tions ( tNavarroetal.il 19961 Il99l hereafter NFW). Fig. [3] shows 
density profiles for a further 9 numerically well-resolved haloes 
from the series of P-WDM runs, with masses ranging from 2.1 x 
lO u h _1 M to 1.5 x lO 12 h _1 M . All nine haloes exhibit well- 
resolved cores and their density profiles are well fit by the pseudo- 
isothermal profile for radii well within the core radius r c out to 
i?200- The isothermal profile fits the halo profiles near R200 to bet- 
ter than 20%, and even better further in. Best-fit values for the cen- 
tral core density vary between the haloes by a factor ~ 3. 

Fig. E] shows the phase space density profiles of the nine P- 
WDM haloes represented in Fig. [3] The blue dashed lines show 
median values in radial bins of the six-dimensional coarse-grained 
phase space density profiles, Fnb, calculated using EnBiD. The 
horizontal red solid lines indicate the theoretically expected max- 
imum fine-grained phase space density, ,/fe> x , from Eqn. (|2j for 
mc 2 = 0.03 keV, the WDM mass appropriate to the initial intrin- 
sic velocities imparted to the particles in the simulations. Clearly, 
all the P-WDM haloes have an approximately flat phase space den- 
sity profile near the centre. The maximum of this coarse-grained 
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Figure 3. Spherically averaged density profiles of nine well-resolved haloes in our fiducial P-WDM simulations (blue triangles). The top three haloes are 
chosen from the higher-resolution P-WDM512 simulation, while the rest are from the three P-WDM256 simulations. The virial mass, M200, in units of 
f0 10 h — 'Mq, is indicated in each panel. Fits using the pseudo-isothermal profiles of Eqn. {5) over the radial range 2e < r < rzoo are shown by long dashed 
lines, with the best-fit core radius, r c , indicated by vertical red lines. For comparison, the best fit NFW profile over the radial range r c < r < T200 is shown 
as the black dot-dashed line. The lower plots show the ratio of the measured densities to the best-fit NFW and pseudo-isothermal profiles. All nine haloes show 
a well-resolved core, and their profiles are well approximated by the pseudo-isothermal form from well inside r c out to -R200 ■ 



density is very close to the fine-grained bound appropriate for the 
mc 2 = 0.03 keV WDM model. 

The red horizontal dashed lines in Fig.[4]show -F ( ™ ax , the max- 
imum coarse-grained phase space density estimated from Eqn. {5} 
assuming a pseudo-isothermal profile and Maxwellian velocities. 
To calculate F;™ ax we used the value of r c which best fits the den- 



sity profiles of Fig[5] and computed the one-dimensional velocity 
dispersion, a, for all particles within r c . For all nine haloes, the 
coarse-grained estimate, Fj™ ax , is within 30 per cent of /™. Thus, 
it is possible to estimate /™ by measuring the core radius and ve- 
locity dispersion of the halo. If real galaxy haloes are dominated by 
WDM and baryons have little effect, then our simulations suggests 
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Figure 4. Spherically averaged coarse-grained phase space density profiles for the nine well-resolved haloes shown in Fig|5] The blue dashed lines show the 
median values in bins in radius, of the six-dimensional coarse-grained phase space density profiles, Fnb, calculated using EnBiD. The pseudo phase space 
density estimate, Q (Eqn ll It . is shown in green, and a power law, F <x r — 19 is shown in black. The red solid line indicates the upper limit to the fine-grained 
phase space density, f^ x , from Eqn. for m<? = 0.03 keV, the value of m adopted in assigning WDM velocities in the simulations. The red dashed 
line indicates the upper limit Fj™ 3 * on F calculated from Eqn. (5) using the value of r c that best fits the density profiles of Fig[3]and the velocity dispersion 
measured for the particles within r c . These values are very close together. The maximum value reached by Fnb as computed using EnBiD (blue dashed line) 
is close to but smaller than the limits calculated by either the coarse-grained estimate Fj™ ax calculated from the halo's properties, or the fine-grained phase 
space limit, /£5J X , appropriate for this P-WDM model. This demonstrates that the simulations satisfy Liouville's theorem, Fnb < /pif x . an d mat the WDM 

particle properties can be estimated from the halo profile, since F™g x ~ ^i™ ax - m contrast, the pseudo phase space density, Q, overestimates Fnb by a 
factor of a few for some haloes, and up to an order of magnitude for others. 
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that it is possible to constrain the mass of a WDM particle, at least 
in principle. We apply this idea to dwarf galaxy data in the next 
section. 

The pseudo phase space density, Q = p/(v 2 ) 3 ^ 2 , often used 
as a proxy for phase space density, is shown with green lines in 
Fig- 111 The shape of Q(r) is similar to that of Fnb(?")» but its am- 
plitude is offset to larger values by up to an order of magnitude. 
Clearly, using Q to estimate _Fnb will result in incorrect constraints 
on WDM particle masses. 



4 CONSTRAINING WDM PARTICLE MASSES FROM 
HALO PROFILES 

We now apply the results of the previous section to the dwarf 
spheroidal satellites of the Milky Way. There are numerous 
claims in the l iterature that t hese g alaxies have a central den- 
sity cor e (e.g. iGilmore et aT] 120071: IWalker & Penarrubial I2012L 
but see IStrigari et alj i 2010bl) ) and there have been attempts 
to interpret these using phase spa ce density argumen t s, of- 
ten with reference to WDM (e.g. Hogan & Dalcanton 2000; 

Bode et alj 1200 ll: IStrigari et al J 12004 



Dalcanton & Hoganl l200ll; 1 

Simon & Gehd 120071: iBovanovskv et alj 120081; IBovarskv et al.l 



2009al : lde Vega& Sanchez 2010; Maccio et al. 20121). These stud 



ies usually assume that the stars in these galaxies are in dynamical 
equilibrium and closely trace the dark matter distribution in their 
parent dark matter haloes. Constraints on the mass of the WDM 
particle are then obtained by comparing the inferred phase space 
density with a theoretical expectation, often b ased on a proxy such 
as th e pseudo phase space density, Q (e.g. iHogan & Dalcanton! 
2000); as we saw before this choice would lead to an incorrect re- 
sult. 

Our analysis shows that the maximum value of the coarse- 
grained phase space density in our N-body simulations is close both 
to the fundamental fine-grained bound, ,/fd x (Eqn. [2}, and to the 
estimate, -F;™ ax (Eqn.|5j, calculated assuming a pseudo-isothermal 
profile and Maxwellian velocities. This result allows us to reassess 
previous constraints on WDM models and to set new rigorous lim- 
its on the WDM particle mass. In particular, we use the results from 
our simulations displayed in Fig. [4] that show that F;™ ax differs 
from /fd x by not more than about 30%. We proceed as follows. 

We assume that the observed half-light radius of a dwarf 
spheroidal is approximately equal to the radius, Rh, at which the 
projected dark matter density attains half its central value. The 
projected surface density profile, S(R), of the pseudo-isothermal 
model can be written as: 



S(R) 



p[(R 2 + z 2 ) 1/2 ]dz 



Sor c 



Vr 2 + R 2 



(17) 



where R is projected radius, and So = npor c is the central sur- 
face density. Therefore R h — V%r c . We estimate the central 
one-dimensional velocity dispersion a of the dark matter from the 
(measured) velocity dispersion of the stars within R^. Combining 
these, we calculate the central maximum phase space density of 
the dark matter halo, F^ x , using Eqn. {3). Finally demanding that 
-FSf ~ /fd x . Eqn. Q, yields an estimate of the WDM particle 
mass, which we expect to be good to within 30 per cent. The in- 
ferred values are given in the last column of Table 2 and plotted in 
Fig. [5] Of course, these constraints are only relevant if the satellite 
dynamics are collisionless, so that Liouville's theorem is satisfied. 

Our constraints on the maximum phase space density and 
mcpg 1 



ferent methods for estimating the maximum value of the coarse- 
grained phase space density. These estimates vary by factors of 
~ 50 to ~ 0.2 for F, depen ding on the method us ed. For example, 
the limits from the model of Hog an & Dalcanton! ( 12000 *) are a fac- 
tor of 3.08 larger than ours. This disagreement is carried over to the 
limits set on mc 2 g 1 ^ A even when using the same data, by factors 
ranging fro m 0.65 to 2.75 (see also the illuminating summary and 
discussion in Bovarsk v et alj d2009al)). In particular, our boun ds on 
m are 1.54 times larger than those of IBovarskv et alj (2009a). Our 
method has the advantage over previous methods that it has been 
explicitly verified by iV-body simulations. 

The method outlined above yields an estimate of the mass of 
the WDM particle, mc 2 (g /2) 1/4 » 0.5 keV, not a limit on m. This 
estimate is based on phase space considerations which, in turn, de- 
pend on the initial intrinsic velocities of the WDM particles. The 
lower limit inferred from the Lyman-a forest, mc 2 ~ 1.2 keV, 
on the other hand, is determined by the cutoff in the power spec - 
trum dSeliaketalJl2006l: Iviel et alj|2008l; IBovarskv et ai] |2009a). 
Although this cutoff is, of course, induced by free streaming due to 
the initial velocities, the two mass estimates exploit different prop- 
erties of the initial WDM particle distribution, velocities in one case 
and the shape of the power spectrum in the other. For two of the 
satellites listed in Table [2]- Coma Berenices I and Canes Venat- 
aci II - the mass estimate from the phase space constraint and the 
lower limit from the Lyman-a forest are close. However, for the 
majority of the dwarf spheroidals in the table, the WDM particle 
mass inferred from the phase space constraint is significantly be- 
low the lower limit from the Lyman-a forest. For these objects the 
two methods can only be reconciled if the phase space density is 
lower than predicted by Liouville's theorem. This may result, for 
example, from energy exchange between the dark matter and the 
baryons in the halo. 



5 SUMMARY AND DISCUSSION 

Whether warm dark matter particles decouple as thermal relics or 
form from non-equilibrium decay, they acquire initial velocities 
whose amplitude depends on the particle's mass. Subsequent free 
streaming imposes a cut-off in the primordial spectrum of density 
perturbations. In this paper we have performed a series of numer- 
ical experiments to investigate how the intrinsic primordial veloc- 
ity dispersion of fermionic dark matter particles affects the central 
density profile of the dark matter haloes into which they later col- 
lect. For WDM the initial velocities are small and resolving them 
in an iV-body simulation of halo formation would require a cur- 
rently prohibitively large number of simulation particles. Since we 
are primarily interested in the connection between the initial veloc- 
ities and the final phase space distribution of the particles, we can 
circumvent this problem by decoupling the initial velocities from 
the free streaming length. Our simulations therefore do not corre- 
spond to a self-consistent representation of any particular WDM 
particle candidate but they are suitable for tackling the problem in 
hand. The power spectrum cut-off and primordial velocities we as- 
sumed correspond to those of thermally produced WDM particles 
of mass of 2 keV and 0.03 keV respectively. For comparison pur- 
poses, we also ran simulations of the standard ACDM case and a 
self-consistent WDM model with a particle mass of 2 keV. 
Our main results may be summarised as follows: 

(i) Initial particle velocities induce cores in the radial profiles of 
both the physical and the phase space density of dark matter haloes. 
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Figure 5. WDM particle mass, mc 2 (g/2) 1 / 4 , inferred from the satellite properties of Table 2, as a function of the satellite velocity dispersion, a (left), and 
core radius, r c (right). For clarity, uncertainties on a and r c are not plotted. The value of m assumes that the coarse-grained phase space density, i^™ ax 
(Eqn.[3}, equals the theoretical bound, /jj 1 ™ (Eqn.fT}. This analysis gives val ues of mc 2 (q/2) 1 / 4 0.5 keV, approximately independently of a or r c . 
Hashed regions indicate the lower bound on m from the Lyman-o forest study of Bovars kv et al . 1 2009a). with their 99.7% and 95% confidence intervals (1.5 
and 1 .7keV) shown in red and green respectively. 



The inner density profile of the simulated haloes is well described 
by the pseudo-isothermal model (with a core). 

(ii) The maximum coarse-grained phase s pace density of simu- 
lated haloes (computed using the EnBid code; ISharma & Steinmetzl 
l2006t) is very close to the theoretical fine-grained upper bound. This 
implies that it is, in principle, possible to use phase space arguments 
to constrain the nature of the dark matter. 

(iii) In contrast, the pseudo phase space density, Q ~ p/o" 3 , 
overestimates the maximum phase space density by a significant 
factor, up to an order of magnitude. 

(iv) Assuming that the velocity distribution of the halo dark mat- 
ter particles is Maxwellian, a simple analytical model that pre- 
dicts the maximum allowed coarse-grained phase space density de- 
scribes the simulations remarkably well (Fig. 3). 

(v) Application of this analytic model to the kinematical data 
of dwarf spheroidal satellites of the Milky Way, assumed to 
have a pseudo-isothermal density profile with a core, constrains 
the mass of a hypothetical thermal fermionic WDM relic to be 
mc 2 (ff/2) 1/4 w 0.5 keV. 

The low particle mass we infer from phase space consider- 
ations yields a large free-str eaming mass, Mfs ~ 10 12 ft _1 M© 
(Eqnll6t. As emphasized by iMaccio et afl j2012l) . this is so large 
that a model with this kind of dark matter is unlikely to produce 
enough satellites in galaxies like the Milky Way. Whereas CDM 
may suffer from an excess of massive subhaloes - the 'too big to 



fail' problem of iBovlan-Kolchin et all fcoilh . WDM may suffer 
from a 'too small to succeed' problem. 

Our inferred value for m differs from other estimates in the lit- 
erature based on different methods to estimate the coarse-grained 
phase space density, by factors ranging from 0.2 to 58. Our esti- 
mate, however, is the only one that has been explicitly validated 
using cosmological simulations of halo formation. 

The inferred value of m follows from the assumption that 
the density profiles of the Milky Way dwarf spheroidal satel- 
lites h ave central cores. Thi s assu mption is controversial: for ex- 
ample Walk er & Pefiarrubial <2012h argue that cores are required 
by the data, at least for Fornax and Sculptor, but Strig ari et al.l 
(2010b) have shown explicitly that the data for these and other 
dwarfs are consistent with NFW cusps. In any case, the value 
of the particle mass required by the kinematical data, under the 
assumption that the dwarf spheroidal haloes do have cores, con- 
flicts with the lower bound on the WDM particle mass set by ob- 
servations of the Lyman-a forest which require the particles to 
have a mass m ~ 1 2 keV dSeliak et alj|2006t IViel etafl l2008t 
iBovarskv et al.ll2009al) . This implies that if the cores are actually 
real, then they cannot be explained by the free-streaming velocities 
of thermally produced WDM particles. Instead, baryonic processes 
associated with th e forming galaxy, for example, of the kind origi- 
nally proposed bv lNavarro et al, Ul996b and more recently seen in 
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Table 2. Parameters for dwarf spheroidal satellites of the Milky Way: satellite name (col. 1); half-light radius and velocity dispersion compiled by 
Boyar sky et alj J2009al) (cols. 2 and 3); coarse-grained limit, ^™ ax , derived from the values of and a (col. 4); mass of the WDM particle assuming 
that the coarse-grained phase space, i^™ ax , equals the fine-grained value, fp^ x (col. 5). 



simulations (e.g.lRead & Gilmordl2005l iMashchenko et alj|2008l : 
iGovemato et al.l20ld . lBrooks & Zolotovl2012h would be required. 
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